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A common  problem  faced  by  simulators  is  that  of  constructing  a 
confidence  interval  for  the  steady-state  mean  of  a stochastic  process. 
We  have  reviewed  the  existing  procedures  for  this  problem  and  found 
that  they  all  produce  confidence  intervals  with  coverages  which  may 
be  considerably  lov;er  than  desired.  Thus,  in  many  cases  simulators 
will  have  more  confidence  in  their  results  than  is  justified. 

In  this  paper  we  present  a new  sequential  procedure  based  on 
the  method  of  batch  means  for  constructing  a confidence  interval 
with  coverage  close  to  the  desired  level.  Empirical  results  for  a 
large  number  of  stochastic  systems  indicate  that  the  new  procedure 
performs  quite  v/ell. 


*This  research  was  supported  by  the  Office  of  Naval  Research  under 
contract  N00014-76-C-0403  (NR  047-145). 


Let  [X.,i  > 1}  be  a stochastic  process  with  steady-state  mean 

t — 


li  = lim  ^ Kjn  (with  probability  1) 


n -*•«> 


i=l 


A common  problem  faced  by  simulators  is  that  of  constructing  a confi- 
dence interval  (c.i.)  for  u-  There  are  two  basic  approaches:  (1) 

Construct  a c.i.  from  an  arbitrary  fixed  sample  size  (2)  Sequentially 
increase  the  sample  size  until  an  "acceptable"  c.i.  can  be  constructed. 

For  the  fixed  sample  size  approach,  five  methods  have  been  sug- 
gested in  the  simulation  literature:  replication,  batch  means,  spectrum 

analysis,  autoregressive  representation,  and  regeneration  cycles.  (See 
Crane  and  Iglehart  [6,7],  Fishman  [10,11],  Iglehart  [12],  and  Law  [14].) 
Unfortunately,  all  of  these  methods  have  the  drawback  that  if  the  total 
sample  size  is  chosen  too  small,  then  the  actual  coverage  of  a con- 
structed c.i.  may  be  considerably  lower  than  desired.  This  was  empir- 
ically shown  for  replication  and  batch  means  in  [14],  and  will  be 
demonstrated  for  the  other  three  methods  in  Law  [16].  (These  results 
v/ere  reported  at  the  1976  Winter  Simulation  Conference.) 

If  the  total  sample  size  is  made  sufficiently  large,  then  the 
actual  coverace  will  be  close  to  the  desired  level;  hov/ever,  what  is 
sufficiently  large  for  one  stochastic  system  may  not  be  adequate  for 
another.  Thus,  in  practice  a simulator  will  not  know  the  actual 
coverage  of  his  c.i . 
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The  above  summary  suggests  that  a procedure  is  needed  to  determine 
the  total  sample  size  necessary  to  achieve  "acceptable"  coverage.  Three 
such  procedures  have  been  suggested.  Fishman  [9]  applied  a procedure 
based  on  autoregressive  representation  to  100  independent  simulation 
runs  of  an  A///V/1  queue  with  p = .9.  On  each  run  he  attempted  to  con- 
struct a 90%  c.i.  for  the  mean  number  of  customers  in  system,  L.  He 
found  that  between  66  and  79  percent  of  the  c.i.'s  covered  L,  depending 
on  the  choice  of  the  initial  sample  size.  Robinson  [20]  applied  a 
similar  sequential  procedure  based  on  regeneration  cycles  to  100  runs 
of  the  W/W/1  queue  with  p = .5.  On  each  run  he  attempted  to  construct 
a 90%  c.i.  for  the  mean  delay  in  queue,  d.  He  found  that  between  60 
and  63  percent  of  the  c.i.’s  covered  d,  depending  on  the  initial  sample 
size.  Mechanic  and  McKay  [17]  developed  a procedure  based  on  batch 
means  that  was  difficult  to  understand  and  was  never  substantially  ,s 

tested.  However,  it  was  their  work  that  provided  the  initial  motivation 
for  the  sequential  procedure  presented  in  the  next  section. 

The  remainder  of  this  paper  is  organized  as  follows.  In  Sections 
1 and  2,  respectively,  we  describe  and  justify  the  new  procedure.  The  ' 

results  of  testing  the  procedure  on  a large  number  of  queueing,  inven- 
tory,  and  computer  models  are  given  in  Section  3.  Finally,  in  Section 

1 

4 we  summarize  our  findings  and  comment  on  the  computational  efficiency 
of  the  procedure. 

1 . The  Procedure 

Suppose  we  make  a simulation  run  of  length  n and  then  divide  the 
resulting  observations  . . . ,X^  into  k batches  of  length  n (n=k-m). 

Let  X.(m)  ( j=]  ,2 , . . . ,k)  be  the  sample  mean  of  the  rr>  observations  in  the 

J 

I 

A 
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jth  batch  and  let  X(k,m)  = ^ Xj(mJ/k  be  the  grand  sample  mean.  If  m 

is  sufficiently  large,  then  the  X.(m)‘s  will  be  essentially  uncorrelated 

J 

(see  Section  2)  and  an  approximate  100(1  - a)%  c.i.  for  p is  given  by 


X(k,m)  ^ \.X(k,m)  j , 


(1) 


where  l-a/2  l-c»/2  point  for  a t distribution  with  k-1  degrees 

of  freedom  and 


(2) 


This  approach  to  constructing  a c.i.  is  called  the  method  of  batch  means 
(see  [14]). 

The  validity  of  the  c.i.  given  by  (1)  depends  crucially  on  the  x.M's 

J 

being  approximately  uncorrelated  (see  Section  2).  We  will  attempt  to 

determine  the  presence  of  significant  correlation  by  estimating  Pjf'm;, 

the  lag  1 correlation  between  the  X.fmJ's.  The  usual  estimator  of  p^{nJ 

«7 

is 

k-1  k 2 

P^(k,m)  = ^ \xym)-f(k,m) 

^ J = 1 

However,  if  p^(k/2,m)  and  p^(kl2,m)  are,  respectively,  the  usual  lag  1 
estimators  based  on  the  first  k/2  and  last  k/2  batches  (k  is  assumed 
to  be  even),  then  we  can  also  eetimate  p^(m)  by  the  jackknifed  estimator 


X^.^^(m)-X(k 


,m)  XAm)-X(k,m)  . 


V 
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p (k,m)  - 23  (k,m)  - [p^(k/2jfn)  + p^('k/2,mj]/2. 
Ill  I 

We  will  use  p^(k,m)  rather  than  p^(k,m)  to  estimate  p^(m)  since,  in 
general,  it  vn’ll  be  less  biased  (see  Miller  [18]). 

We  now  state  our  sequential  procedure. 

Step  0.  Let  a - .225,  = 600,  = 800,  and  i = 1. 

Step  l.a.  Divide  the  n.  observations  into  400  batches  of  size 

V 

m = n./400.  Compute  p,(400,m)  from  X .(m)  (j=l ,2, . . . ,400) . 

If  pj(400,m)  ^6*,  go  to  Step  l.c.  If  Pj(400,m)  ;<  0 (see 
Note  1),  go  to  Step  2.  Otherwise,  go  to  Step  l.b. 

b.  Divide  into  200  batches  of  size  2m.  Compute  Pi(200,2m) 
from  I. (2m)  (J=l ,2, . . . ,200) . (See  Note  2.)  If  p,(200,2m)  < 
Pj(400,m)  (see  Note  3),  go  to  Step  2.  Otherwise,  go  to 
Step  l.c. 

c.  Replace  t by  t + 1 , set  n.  = 2n.  „ (see  Note  4),  collect 
the  additional  observations  required,  and  go  to  Step  l.a. 

Step  2.  Divide  n.  into  40  batches  of  size  10m.  Use  7. (10m) 

t 0 

(j=l  ,2, . . . ,40)  (see  Note  2)  in  (1)  to  construct  a c.i. 
for  p. 

Notes: 

1.  If  p^(m)  £ 0,  then,  as  is  discussed  in  Section  2,  40  batches  of 
size  10m  will  most  likely  produce  a c.i.  with  at  least  100(1  - a)% 


coverage. 
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2.  An  appropriate  number  of  the  X.{m)'s  may  be  averaged  to  compute 

the  7.(2m)'s  or  the  7.{10m)'s. 

J 0 

3.  If  0 < pj(m)  < .225  and  if  is  decreasing  for  m'  then 

batches  of  size  10m  are  approximately  uncorrelated  (see  Section  2). 
We  test  to  see  if  Pj(m')  is  decreasing  for  m'  by  checking  to 
see  whether  Pj(2m)  < pj(m). 

4.  The  successive  sample  sizes  considered  are  800,1200,1600,2400,...  . 
Thus,  the  total  sample  size  is  doubled  every  other  iteration. 

The  sequential  procedure  which  is  described  above  uses  400  batches 
of  size  m to  decide  when  40  batches  of  size  lOn  are  uncorrelated.  By 
using  400  rather  than  40  batches,  the  correlation  estimator,  which  is 
used  to  determine  the  stopping  point,  has  a smaller  bias  and  variance. 

In  fact,  the  procedure  will  not  work  well  if  the  correlation  is  estimated 
directly  from  40  batches. 

The  next  section  gives  a justification  of  the  procedure.  The  reader 
who  is  primarily  interested  in  the  procedure's  performance  may  proceed 
directly  to  Section  3. 


2.  Justification 

In  the  following  two  subsections  we  discuss  the  general  form  of  the 
sequential  procedure  and  the  choice  of  the  stopping  value  q. 

A.  General  Form  of  the  Procedure 

Suppose  that  the  observations  are  from  a covariance 

stationary  process.  For  t=0,l , . . . ,n-l , let  C.  = Cov[x.tX.  .1  and  for 

^ j j + t 

t»0,l , . . . ,/c-l , let  C.(m)  = Cov[X  .(m)  ,X . .(m)']  and  p.(m)  = C .(m)/C  (m) . 

to  a7  ^ to  0 


The  following  lemma,  which  is  proved  in  the  appendix,  shows  that 


6 


pjm)  : 0 (t  ^ 1 ) for  sufficiently  large  m. 


Lemma  1 . If  0 < then  p.(m)  -*  0 as  m for  i = ] ,Z , . . . ,k-] 


Let  bCk.rn)  be  defined  by 

E[o^[X(k:,m)^)  - bik,m]o'^[X(ktn})'\, 

where  a^\l(k,m)'\  was  given  by  (2).  We  now  prove  that  a^{x(kjm)'\  is 
asymptotically  unbiased  as  m *- 


Theorem  2.  If  0 ^ i;  then  L'k,m)  1 as  m + ” for  all  /c  _>  2. 

1-  - CO 


Proof:  It  is  easy  to  show  that  (see  [14J) 


b(k,m)  = 


k - 1 


The  desired  result  follows  since  -►  0 as  m > for  i=l  ,2 . • • • .*:-l 

by  the  lemma. 

The  above  expression  for  b(k,m)  shows  that  "^^\t(k,m)'\  is  unbiased 
when  QAm)  = 0 for  £ = 1 ,2.  • • • .^'1 . has  a negati'^e  bias  [hik,m)  ■ 1)  when 
p.(m)  > OVi,  and  has  a positive  bias  {b(k,n>  > 1)  when  p.(m>  0 V:. 

t i 

The  case  p^(rrj  > 0 )/ L is  of  the  greatest  concern  since  will 

then  be  underestimated  and  the  coverage  of  the  resulting  c.i.  is  likely 
to  be  less  than  desired. 
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If  we  choose  m large  enough  so  that  the  X.(m)'s  are  approximately 

0 

normally  distributed  in  addition  to  being  uncorrelated,  then  it  becomes 
plausible  to  proceed  as  if  the  X.(m)'s  viere  independent  identically 

t7 

distributed  normal  random  variables  (r.v.'s)  and  to  use  (1)  to  construct 
a c.i . for  p. 

There  are  three  potential  sources  of  error  when  using  (1)  to  con- 
struct a c.i.: 

(1)  Bias  in  o^\_x(k,m)'\  when  m is  too  small  for  the  X.(m)'s  to  be 

J 

uncorrelated. 

(2)  Nonnormality  of  the  XAm)'s. 

(3)  The  fact  that  [X.,i  >^1}  is  not,  in  practice,  covariance  sta- 
tionary. 

However,  for  simple  queueing  models  (e.g.,  Law  [14]  found  that 

the  bias  in  a'^\,X(k,m)']  was  the  most  serious  source  of  error  and  that 
nonnormal ity  was  not  a problem  for  k approximately  20  or  more.  This 
suggests  that  a sequential  procedure  based  on  batch  means  must  be 
able  to  determine  that  batch  size,  m,  for  which  the  X.(m)'%  are 

d 

approximately  uncorrelated. 

To  determine  the  types  of  correlation  which  can  occur  in  practice, 
we  studied  the  following  processes  for  which  Q.(m)  and  h(k,m)  can  be 
analytically  computed: 

(1)  {D.,i  > 1}  for  the  M/M/\  queue  (see  Daley  [8])  with 

p = .5, .8,  and  .9,  where  D.  is  the  delay  in  queue  of  the  tth 

t 

customer. 

« 

(2)  [E.,i  > 1}  for  an  (s,S)  inventory  system  (see  [14)  for  details), 

1 ~ 

where  is  the  expenditure  in  the  tth  period. 


(3)  Thirty  different  AR(1),  AR(2),  and  ARMA(l.l)  time  series 
models  (see  Box  and  Jenkins  [3,  p.  46])  with  parameters 
chosen  over  the  entire  range  of  feasible  values. 

From  studying  these  34  stochastic  processes,  we  found  essentially 
three  types  of  behavior  for  pi(m)  as  a function  of  examples  of 
which  are  shov;n  in  Figures  1,  2,  and  3.  For  type  1 behavior  the 
lag  1 correlation  p^(rn)  strictly  decreases  to  zero.  If  for  some  m, 
PjfmJ  < .4,  then  .9  < ^(40,10^)  < 1 and  pj(lOm)  ~ .05.  That  is,  if 
pjfmJ  < .4,  then  the  variance  estimator  based  on  40  batches  of  size 
10m  is  approximately  unbiased.  The  /■/////!  queue  exhibits  type  1 
behavior. 

In  type  2 behavior,  p^(m)  changes  direction  one  or  more  times 
and  then  strictly  decreases  to  zero.  If  for  some  m,  p^(m)  < .4  and 
p^Cm'J  is  decreasing  for  m'  ^ m,  then  .9  < i>(40,10n)  < 1 and 
Pj(lOm)  = .05.  The  M/M/1  queue  with  service  in  random  order  (SIRO) 
is  of  this  type  (see  Figure  4). 

For  type  3 behavior,  PiCrriJ  < 0 and  2) (40,1  CW;)  > 1,  for  all  m.  In 

this  case  the  1 .nOmJ 's  may  be  correlated,  but  the  coverage  will  be  at 
J 

least  as  great  as  that  desired.  The  (s,S)  inventory  system  exhibits 
type  3 behavior. 

We  certainly  do  not  claim  that  the  above  three  types  of  behavior 
are  the  only  ones  that  can  occur.  In  fact,  for  some  of  the  time  series 
models  we  studied,  p^fmJ  can  be  positive  or  negative.  Hov/ever,  if 
PjfmJ  > 0 for  some  m,  then  we  found  that  type  1 behavior  is  followed, 
and  if  p^(m)  < 0 for  some  m,  then  behavior  similar  to  that  of  type  3 
is  followed. 
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FIGURE  CAPTIONS 


Figure  1.  Type  1 Behavior. 


Figure  2.  Type  2 Behavior. 


Figure  3.  Type  3 Behavior. 


Types  1 and  2 behavior  described  above  explain  why  the  procedure 
checks  for  p^(m)  < a and  p^(m')  decreasing  for  m'  Type  3 behavior 

explains  why  it  checks  for  p^(m)  £ 0. 


B.  The  Choice  of  the  Stopping  Value  a 

i 

% 

! The  above  discussion  was  predicated  upon  knowing  p^(m)  which,  in 

h fact,  is  estimated  by  PjC400,m;.  Thus,  it  is  possible  for  p^(m)  to  be 

' I 

much  larger  than  .4,  but  to  have  the  estimate  Pj('400,m;  < .4,  which 

might  result  in  the  procedure's  stopping  prematurely.  To  determine  how 
i 

small  £7  should  be  to  account  for  the  sampling  variability  of  p^(W0,m), 
we  applied  our  sequential  procedure  with  various  values  of  c to  the 
M/M/}  FIFO  queue,  the  m/m/}  LIFO  queue,  and  the  M/M/}  SIRO  queue,  each 
with  p = .8.  For  each  system  we  made  200  independent  simulation  runs  of 
the  process  {D.,i>}}  and  attempted  to  construct  90%  c.i.'s  for  d = 3.2. 

t — 

The  results  of  these  simulations  are  given  in  Table  I.  Note  that  for 
the  FIFO  and  SIRO  disciplines,  a = .375  gives  good  coverage,  but  that 
for  the  LIFO  case  a much  smaller  value  of  a is  required.  We  chose 
a * .225  as  our  stopping  value  because  we  felt  that  a true  coverage 
of  80%  is  the  minimum  acceptable  coverage  for  a "90%  c.i."  A smaller 
value  of  c would  result  in  better  coverage  for  LIFO,  but  also  in  a 
considerable  increase  in  sample  size  for  the  other  models.  (These 
additional  data  would  not  be  wasted  if  a smaller  c.i.  half  length  was 
desired. ) 

The  fact  that  a smaller  value  of  c?  is  required  for  the  LIFO  model 
is  due  not  to  the  variability  of  pj('400,mJ,  but  rather  to  PjfmVs  being 
a more  slowly  decreasing  functipn  of  m than  for  any  of  the  other  44 
models  considered  in  this  paper.  This  can  be  seen  in  Figure  4 where 
we  plot  p^(m)  (the  average  of  the  200  values  of  p,(400,m))  as  a function 
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TABLE  I.  Proportion  of  Coverage  and  Average 
Sample  Size  for  Various  Values  of  a. 


M/M/1 
P = 

FIFO 

.8 

M/M/1 

P = 

LIFO 

.8 

1 

M/M/1 
P = 

SIRO 

.8 

1 

Q 

Coverage 

Average 

Sample 

Size 

Coverage 

Average 

Sample 

Size 

Coverage 

Averagf 

Sample 

Size 

.375 

.900 

36464 

.650 

3882 

.880 

35616 

.250 

.900 

59776 

.760 

15054 

.860 

59680 

.225 

.890 

67584 

.805 

19100 

.865 

69280 

.200 

.890 

1 

77952 

.825 

26796 

.860 

77472 

.150 

.910 

100992 

.840 

41058 

.865 

99008 

FIGURE  CAPTION 


Figure  4.  p^(m)  as  a Function  of  m for  the  M/M/] 

FIFO,  LIFO,  and  SIRO  Queues,  each  with 

p = .8. 


. Lif-"- 


trtrts  lower  than  FIFO  and 


ot  m for  the  tnree  jjeeei.  Observe 

SIRO  but  decreases  more  slowly.  In  particular,  p^(mJ  < .4  for  m = 6, 
but  Pj(l(>n)  s .19  iPi(m)  < .4  implies  Pi(lOm)  « .05  for  the  other  models). 
By  choosing  a = .225,  we  are  designing  the  sequential  procedure  to  per- 
form adequately  in  the  "worst"  case. 

3.  Empirical  Results 

In  order  to  see  how  well  the  sequential  procedure  works,  we  simulated 
a large  number  of  stochastic  systems  for  which  analytical  results  are 
available.  The  results  of  these  simulations  are  presented  in  the  next 
three  subsections. 

The  random  numbers,  {u.^i  ^1),  used  in  these  simulations  were 
generated  from  the  following  generator  which  is  available  on  the 
Univac  1110: 


J.  ■ (5**r.  , + 1)  mod  2”  (i=l,2,...) 

% I 

u.  » ir./2”  (i=l,2,...), 

% % 

where  y,  is  a specified  seed.  For  a discussion  of  this  generator,  see 
Coveyou  and  Macpherson  [5]. 

A.  Queueing  Systems 

We  first  considered  a variety  of  simple  queueing  systems.  Let 
M,  and  H2  denote,  respectively,  the  4-Erlang,  the  exponential,  and  the 
hyperexponential  distributions  with  coefficients  of  variation  .5,  1, 
and  2.  (See  Law  [13]  for  further  discussion  of  For  each  of  the 

queueing  systems  M/M/^  (p  * .5),  ff^/W/l  (p  » .8),  (p  = .8), 

M/M/Z  (p  » .8),  (p  = .8),  we  made  100  independent  simulation 

runs  of  the  stochastic  process  {D.,i  > 1}  and  attempted  to  construct 
90*  c.i.'s  for  d.  For  each  system,  E(A)  * 1 (the  mean  interarrival 


14 


time)  and  = 0 (i.e.,  no  customers  are  present  at  time  zero).  In 
Table  II  we  give  for  each  system  d,  the  proportion  of  the  100  c.i.'s 
which  covered  d,  the  average  sample  size  at  termination,  and  the  aver- 
age c.i.  half  length  divided  by  d.  (We  will  henceforth  call  the  aver- 
age c.i.  half  length  divided  by  p the  "relative  precision"  of  the  c.i.) 
For  completeness  we  also  give  the  previous  results  for  the  A//W/1  queue 
with  p = .8  and  FIFO,  LIFO,  and  SIRO  queue  disciplines.  Observe  that 
the  average  sample  size  at  termination  is  quite  system  dependent.  For 
example,  the  M/W/1  queue  with  p = .8  requires  a sample  size  about  8 
times  larger  than  that  required  by  the  W/W/1  queue  with  p = .5.  This 
is,  of  course,  due  to  the  fact  that  the  process  {D^,i  1}  is  more  cor- 

related for  larger  values  of  p. 

For  the  W/W/1  queue  it  can  be  shown  that  (see  [8])  an  approximate 
theoretical  sample  size  of 


„ . iE(A)yoH2  ^ 5o  - 4p"  p^ 

(1  - p)^ 


(3) 


(1  - P)' 

is  required  to  obtain  a c.i.  whose  half  length  is  6d  (0  < 6 < 1). 

Since  the  average  half  length  for  p = .5  was  .098  of  d (see  Table  II), 
if  we  substitute  6 = .098  into  (3),  we  get  a theoretical  sample  size  of 
8572.  Similarly,  substitution  of  6 = .068  into  (3)  gives  a theoretical 
sample  size  of  75813  for  p * .8.  These  sample  sizes  compare  closely  to 
those  actually  obtained  by  the  procedure. 


B.  An  Inventory  System 

The  second  type  of  stochastic  system  we  considered  was  an  (s,S) 
inventory  system  with  zero  delivery  lag  and  backlogging.  This  system 


TABLE  II.  Empirical  Results  for  Queueing  Systems 


System 

fl 

d 

Number 

of 

Runs 

Proportion 

of 

Coverage 

Average 

Sample 

Size 

Relative 

Precision 

W/W/1  FIFO 

.5 

.50 

100 

.850 

8352 

.098 

W/W/1  FIFO 

.8 

3.20 

200 

.890 

67584 

.068 

W/M/1  LIFO 

.8 

3.20 

200 

.805 

19100 

.141 

MIMI\  SIRO 

.8 

3.20 

200 

.865 

69280 

.067 

.8 

1.81 

100 

.830 

42240 

.075 

MiHzn 

.8 

8.00 

100 

.800 

188928 

.072 

Ml  Ml  2 

.8 

2.84 

100 

.910 

65792 

.076 

W/W/l/W/1 

.8 



6.40 

100 

.900 

88832 

.049 
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is  described  in  detail  in  [14].  We  made  100  independent  simulation 
runs  of  the  process  ^ 1}  and  attempted  to  construct  90%  c.i.'s 

for  the  steady-state  mean  expenditure  per  period,  e = 112.108.  We 
found  that  each  run  terminated  on  the  first  iteration  {n  = 800)  and 
that  99  out  of  100  c.i.'s  covered  e.  These  results  can  be  explained 
by  the  fact  that  p.(m)  *-  0 (i  ^ 1)  ^nd  b(40,10m)  > 1 for  this  system 

(see  [14]  for  some  actual  values). 

C.  Computer  Models 

To  see  how  the  sequential  procedure  works  on  more  complex  systems, 
we  simulated  two  models  of  computer  systems.  These  models  resemble  more 
closely  the  type  of  system  that  might  actually  be  simulated  in  the  real 
world  than  do  the  simple  models  of  the  previous  two  subsections.  We 
first  considered  a model  of  a time-shared  computer  system  which  was 
studied  by  Adiri  and  Avi-Itzhak  [1]  and  more  recently  in  a simulation 
context  by  Sargent  [19].  There  are  N terminals  and  a single  central  pro- 
cessing unit  (CPU)  as  shown  in  Figure  5.  The  operator  of  each  terminal 
"thinks"  for  an  amount  of  time  which  is  an  exponential  r.v.  with  rate 
Pi  and  then  sends  a message  with  service  time  which  is  an  exponential 
r.v.  with  rate  P2-  The  arriving  jobs  join  a FIFO  queue  at  the  CPU. 

The  CPU  allocates  to  each  job  a maximum  service  quantum  of  length  q. 

If  the  (remaining)  service  time  of  a job,  s,  is  less  than  or  equal  to 
q,  then  the  CPU  spends  time  s plus  t (a  fixed  overhead)  processing 
the  job  and  the  job  returns  to  its  terminal.  If  s > q,  then  the  CPU 
spends  time  q plus  t processing  the  job  and  the  job  joins  the  end  of 
the  queue.  This  process  is  repeated  until  the  job's  service  time  is 
eventually  completed. 


HGURE  CAPTIONS 


Figure  5.  Time-Shared  Computer  Model. 


Figure  6.  Central -Server  Computer  Model. 
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We  let  N = 35,  Pj  = 1/25,  = 5/4,  <7  = .1,  t = .015,  and  then 

simulated  the  stochastic  process  [R.,i  ^ 1},  where  R.  is  the  response 
time  of  the  ith  job  requesting  service  and  all  terminals  are  in  the 
think  state  at  time  zero.  Our  objective  was  to  construct  a 90%  c.i. 
for  the  steady-state  mean  response  time,  v - 8.246.  lie  made  100  inde- 
pendent runs  and  obtained  a coverage  of  .92,  an  average  sample  size  of 
36160,  and  a relative  precision  of  .036. 

We  next  considered  the  central -server  model  of  a mul ti programmed 
computer  system  (see  Buzen  [4]).  There  is  a CPU  (unit  1)  and  W-1 
peripheral  units  (units  2 through  M)  as  shown  in  Figure  6.  Each  unit 
has  its  own  FIFO  queue  and  the  service  time  at  unit  i is  an  exponential 
r.v.  with  rate  p^  (T=l  ,2, . . , . It  is  assumed  that  there  are  R jobs  in 

the  system  at  all  times.  When  a job  completes  service  at  the  CPU  it 
leaves  the  system  with  probability  Pj  or  it  goes  to  peripheral  unit  i 
with  probability  p^  (t=2,3,  . . ,A/) . ^ job  leaving  the  system  is  instanta- 
neously replaced  by  a new  job  which  joins  the  end  of  the  CPU  queue.  A 
job  leaving  a peripheral  unit  also  joins  the  CPU  queue. 

We  made  100  simulation  runs  of  the  process  iR^^i  1 D 
different  sets  of  parameters  and  on  each  run  we  attempted  to  construct 
a 90%  c.i.  for  r (the  steady-state  mean  time  between  entries  to  the  CPU 
queue).  The  parameters  for  these  models  and  the  results  of  the  simula- 
tions are  given  in  Table  III.  Also  included  are  the  steady-state  proba- 
bility that  unit  i is  busy,  p^,  and  the  state  of  the  system  at  time 
zero,  - ■ ■ ,l„)  ^ where  1.  is  the  number  of  jobs  at  unit  i. 

M t 
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We  have  used  the  sequential  procedure  to  construct  "90%  c.i.'s" 
for  13  stochastic  systems  (excluding  the  inventory  model),  obtaining 
coverages  between  .80  and  .92  and  relative  precisions  between  .027 
and  .141.  Averaging  over  the  13  systems,  we  obtain  a mean  coverage 
of  .867  and  a mean  relative  precision  of  .067. 

The  sequential  procedure  described  here  is  fairly  easy  to  program, 
is  computationally  efficient,  and  requires  only  800  storage  locations 
to  store  the  batch  means.  A FORTRAN  program  for  the  procedure  and  an 
explanation  of  how  to  use  it  may  be  found  in  Law  and  Carson  [15]. 


Appendix 

Our  objective  is  to  prove  Lemma  1 (see  Section  2), 


m-^ 


Lemma  A.l.  If  ^ converges,  then  ^ 0 as  m 


)l=l 


£=1 


m-1  s m-1  m-1 

Proof;  Y V”  ' Z ■ Z 


(m  - 


£=1 


£=1  £=1 


The  result  follows  since  the  second  sum  on  the  right-hand  side  converges 


to  ^ by  Lemma  8.3.1  of  Anderson  [2,  p.  460]. 


£=1 


Letima  A. 2.  If  ^ converges,  then  SJmJ 
£=1 


m-1 


= Z " - w 


.7  = 1 


as  m -*  «>  for  t > 1 . 
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Proof:  ^ ~ Y V"'  " Z " " Z 


J,=tr77+1 


a=l 


^=1 


l=im+^ 


The  first  and  fourth  sums  converge  to  0 since  I converges.  The  second 


and  third  sums  converge  to  0 by  Lemma  A.l 


I 


Lemma  A. 3.  If  } ,C.  converges,  then  s‘'.(m)  = ) (l-j7m)C.  > 0 as 


£=1 


m~^ 

1 


^=l 

m -*■  oo  for  i > 1 . 


.7  = 1 


m-1  m-1 

Proof:  S^.(m)  = y {im-j)C.  ./m  - (i  - 1)  / C. 

i L-i  tm-j 

J=1  J=1 


im-1 


im-] 


= ^ IC^/m  - ^ IC^/m  - (t  - 1)  Y. 


S- 


il=l 


£=1 


)m+l 


The  first  and  second  sums  converge  to  0 by  Lemma  A.l.  The  third  sum 

00 

converges  to  0 since  ^ c.  converges. 


£=1 


Lemma  1 . I f 0 ■ ^ 

£=  - oo 


i’  ■'  w,  then  p.f'nj  - 0 as 
Z I 


for  : = 1,2 A:-l 


Proof:  Since  P^(m)  = mC  Am) jmC  and  mC^(m)  ^ ^ as  m 


i=  -oo 


idMH 
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’I 


(see  [2,  p.  459]),  it  suffices  to  show  that  mCAm)  0 as  m -*■  ®.  However, 

mC.(m)  = C.  + s\(m)  + S^.(m)  (see  [14])  and  C.  -*■  0 as  m -*■  «>  since 
^ tm  t V ' im 

OO 

^ converges.  The  desired  result  now  follows  from  Lemmas  A. 2 and  A. 3. 
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